# ifndef _ADVREACDIFF_HPP
# define _ADVREACDIFF_HPP 1

/** include predefined life command line options */
#include <life/options.hpp>

/** include linear algebra backend */
#include <life/lifealg/backend.hpp>

/** include function space class */
#include <life/lifediscr/functionspace.hpp>

/** include helper function to define \f$P_0\f$ functions associated with regions  */
#include <life/lifediscr/region.hpp>

/** include integration methods */
#include <life/lifepoly/im.hpp>

/** include gmsh mesh importer */
#include <life/lifefilters/gmsh.hpp>

/** include exporter factory class */
#include <life/lifefilters/exporter.hpp>

/** include gmsh generator for tensorized domains */
#include <life/lifefilters/gmshtensorizeddomain.hpp>

/** include  polynomialset header */
#include <life/lifepoly/polynomialset.hpp>

/** include  the header for the variational formulation language (vf) aka FEEL++ */
#include <life/lifevf/vf.hpp>

#include <life/lifemesh/meshmover.hpp>

/** use Life namespace */
using namespace Life;
using namespace Life::vf;

/**
 * This routine returns the list of options using the
 * boost::program_options library. The data returned is typically used
 * as an argument of a Life::Application subclass.
 *
 * \return the list of options
 */
inline
po::options_description
makeOptions()
{
    po::options_description laplacianoptions("Laplacian options");
    laplacianoptions.add_options()
        ("hsize", po::value<double>()->default_value( 0.5 ), "mesh size")
        ("epsilon", po::value<double>()->default_value( 0.01 ), "mesh size")
        ("betaX", po::value<double>()->default_value( 1 ), "grad.grad coefficient")
        ("betaY", po::value<double>()->default_value( 1 ), "grad.grad coefficient")      
        ("mu", po::value<double>()->default_value( 1 ), "grad.grad coefficient")
        ("stab", po::value<int>()->default_value( 1 ), "ajout des termes de stabilisation" )
        ("weakdir", po::value<int>()->default_value( 1 ), "use weak Dirichlet condition" )
        ("penaldir", Life::po::value<double>()->default_value( 10 ),
         "penalisation parameter for the weak boundary Dirichlet formulation")
        ;
    return laplacianoptions.add( Life::life_options() );
}

/**
 * This routine defines some information about the application like
 * authors, version, or name of the application. The data returned is
 * typically used as an argument of a Life::Application subclass.
 *
 * \return some data about the application.
 */
inline
AboutData
makeAbout()
{
    AboutData about( "MonALE" ,
                     "MonALE" ,
                     "0.2",
                     "nD(n=1,2,3) Laplacian on simplices or simplex products",
                     Life::AboutData::License_GPL,
                     "Copyright (c) 2008-2009 Universite Joseph Fourier");

    about.addAuthor("Vincent Chabannes", "developer", "vincent.chabannes@imag.fr", "");
    return about;

}


/**
 * \class AdvReacDiff
 *
 * Résoud un problème d'advection-diffusion-réaction
 * avec possibilitée de stabilisation
 *
 */
template<int Dim,int Order>
class AdvReacDiff
    :
    public Application
{
    typedef Application super;
public:

    //! Polynomial order \f$P_2\f$
//    static const uint16_type Order = 5;

    //! numerical type is double
    typedef double value_type;

    //! linear algebra backend factory
    typedef Backend<value_type> backend_type;
    //! linear algebra backend factory shared_ptr<> type
    typedef boost::shared_ptr<backend_type> backend_ptrtype;


    //! sparse matrix type associated with backend
    typedef typename backend_type::sparse_matrix_type sparse_matrix_type;
    //! sparse matrix type associated with backend (shared_ptr<> type)
    typedef typename backend_type::sparse_matrix_ptrtype sparse_matrix_ptrtype;
    //! vector type associated with backend
    typedef typename backend_type::vector_type vector_type;
    //! vector type associated with backend (shared_ptr<> type)
    typedef typename backend_type::vector_ptrtype vector_ptrtype;

    //! geometry entities type composing the mesh, here Simplex in Dimension Dim of Order 1
    typedef Simplex<Dim> convex_type;
    //! mesh type
    typedef Mesh<convex_type> mesh_type;
    //! mesh shared_ptr<> type
    typedef boost::shared_ptr<mesh_type> mesh_ptrtype;


    typedef Mesh<Simplex<1, 1> > struct_mesh_type;
    typedef bases<Lagrange</*N*/1, Scalar, PointSetFekete> > struct_basis_type;
    typedef FunctionSpace<struct_mesh_type, struct_basis_type, double> struct_functionspace_type;
    typedef boost::shared_ptr<struct_functionspace_type> struct_functionspace_ptrtype;
    typedef typename struct_functionspace_type::element_type struct_element_type;


    //! function space that holds piecewise constant (\f$P_0\f$) functions (e.g. to store material properties or partitioning
    typedef FunctionSpace<mesh_type, bases<Lagrange<0,Scalar> >, Discontinuous > p0_space_type;
    //! an element type of the \f$P_0\f$ discontinuous function space
    typedef typename p0_space_type::element_type p0_element_type;

    //! the basis type of our approximation space
    typedef bases<Lagrange<Order,Scalar> > basis_type;

    //! the approximation function space type
    typedef FunctionSpace<mesh_type, basis_type> space_type;
    //! the approximation function space type (shared_ptr<> type)
    typedef boost::shared_ptr<space_type> space_ptrtype;
    //! an element type of the approximation function space
    typedef typename space_type::element_type element_type;

    //! the exporter factory type
    typedef Exporter<mesh_type> export_type;
    //! the exporter factory (shared_ptr<> type)
    typedef boost::shared_ptr<export_type> export_ptrtype;




    typedef bases<Lagrange<1, Vectorial, PointSetFekete> > dep_basis_type;
    typedef FunctionSpace< mesh_type, dep_basis_type, double> dep_functionspace_type;
    typedef boost::shared_ptr<dep_functionspace_type> dep_functionspace_ptrtype;
    typedef typename dep_functionspace_type::element_type dep_element_type;

    /**
     * Constructor
     */
    AdvReacDiff( int argc, char** argv, AboutData const& ad, po::options_description const& od )
        :
        super( argc, argv, ad, od ),
        M_backend( backend_type::build( this->vm() ) ),
        meshSize( this->vm()["hsize"].template as<double>() ),
        exporter( Exporter<mesh_type>::New( this->vm(), this->about().appName() ) )
    {
    }

    /**
     * create the mesh using mesh size \c meshSize
     *
     * \param meshSize the mesh characteristic size
     * \return a mesh of an hyper cube of dimension Dim
     */
    mesh_ptrtype createMesh( double meshSize );

    /**
     * run the convergence test
     */
    void run();

private:

    /**
     * solve the system \f$D u = F\f$
     *
     * \param in D sparse matrix
     * \param inout u solution of the system
     * \param in F vector representing the right hand side of the system
     */
    void solve( sparse_matrix_ptrtype& D, element_type& u, vector_ptrtype& F );


    /**
     * export results to ensight format (enabled by  --export cmd line options)
     *
     * \param u the function to save in ensight format
     */
    void exportResults( element_type& u ,element_type& uexact );

private:

    //! linear algebra backend
    backend_ptrtype M_backend;

    //! mesh characteristic size
    double meshSize;

    //! exporter factory
    export_ptrtype exporter;
}; // Laplacian

# endif

